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We study control of synchronization in weakly coupled oscillator networks by using a phase re¬ 
duction approach. Starting from a general class of limit cycle oscillators we derive a phase model, 
which shows that delayed feedback control changes effective coupling strengths and effective fre¬ 
quencies. We derive the analytical condition for critical control gain, where the phase dynamics 
of the oscillator becomes extremely sensitive to any perturbations. As a result the network can 
attain phase synchronization even if the natural interoscillatory couplings are small. In addition, we 
demonstrate that delayed feedback control can disrupt the coherent phase dynamic in synchronized 
networks. The validity of our results is illustrated on networks of diffusively coupled Stuart-Landau 
and FitzHugh-Nagumo models. 


I. INTRODUCTION 

Starting from C. Huygens’ research on “an odd kind 
sympathy” between coupled pendulum clocks, the syn¬ 
chronization as a phenomenon occurs in various man¬ 
made and natural systems 0-0. The coherent behavior 
of oscillators arises in numerous situations, e.g., flash¬ 
ing of fireflies [§}, cardiac pacemaker cells Q neurons 
in the brain [7|, coupled Josephson junctions |8(, chem¬ 
ical reactions |l|, 0, crowd synchrony 0, and power 
grids 0,0. The synchronous behavior can be desirable 
or harmful. The ability to control synchrony in oscilla¬ 
tory networks covers a wide range of real-world applica¬ 
tions, starting from neurological treatment of Parkinsons 
disease and essential tremor 00 to the design of ro¬ 
bust power grids 0,0 . 

Phase reduction is a fundamental theoretical technique 
to investigate sjmchronization in weakly coupled oscil¬ 
lator networks [f]-[4j, since it allows the approximation 
of high-dimensional dynamics of oscillators with a single 
phase variable. The concept of the phase model causes 
significant progress in understanding the synchrony of the 
networks, e.g., correlation between topology and dynam¬ 
ics towards synchronization [16|. synchronization crite¬ 
rion for almost any network top olog y [12j, optimal syn¬ 
chronization 0 , chimera states [Hi Il9j|. etc. The main 
factors determining the synchrony in the phase model are 
coupling strength and dissimilarity of frequencies. The 
ability to change these parameters will easily allow the 
synchronization or desynchronization of networks. Typ¬ 
ically, the phase variable is not attained for direct mea¬ 
surements and actions. Instead of this, we have an access 
to dynamical variables of the limit cycle. In such situ¬ 
ations, the control schemes are usually based on feed¬ 
back loops. Therefore, we ask: how do we enhance or 
suppress synchronization in networks via feedback sig¬ 
nals, when minimal knowledge about the particular unit 
of the network is available? This question has been in¬ 
vestigated in Refs. 0,0], where the synchronization is 
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controlled by delayed mean-field feedback into the net¬ 
work. Also, there have been many numerical investiga¬ 
tions 0-0 devoted to this question. In this work, we 
present analytical results for control of synchronization 
in networks by time delayed local signals fed back into 
particular units of the network. By using phase reduction 
for systems with time-delay [2g, 0], we arrive at a phase 
model. It fully coincides with the phase model of the 
uncontrolled network, the only difference being that the 
coupling strengths and frequencies depend on the con¬ 
trol parameters. Surprisingly, the relations are almost 
universal, i.e. do not depend on the particular model of 
the limit cycle and on the coupling interoscillatory func¬ 
tion. Moreover, the coupling strengths have a multiplica¬ 
tive inverse dependence on the feedback control gain and 
their values can be selected from zero to infinity. As a 
consequence, the synchronization can be achieved even if 
the oscillators are almost uncoupled. Also, we show that 
the particular choice of the delay times in the control 
scheme can lead to full phase synchronization (i.e., when 
the phases of all oscillators are equal at any time mo¬ 
ment) in the network. The analytical results are verified 
numerically on networks of diffusively coupled Stuart- 
Landau and FitzHugh-Nagumo models. 


II. PHASE REDUCTION OF OSCILLATOR 
NETWORK 

We consider a general class of N weakly coupled limit 
cycle oscillators under delayed feedback control (DFC) 

N 

Xj = f,; (x,) + e X G U (*U x l) + K ‘ [ X i(* - Ti) ~ x i W] , 

3 =1 

(1) 

where x^ £ R” is an n-dimensional state vector of the 
ith oscillator, : R" —>• R” is a vector field repre¬ 
senting the free dynamics of the Zth oscillator, G ij : 
R” x R" —>• R” is an inter-oscillatory coupling function, 
Ki = diagf-ff^, K- 2 \ ..., K^] is an n-dimensional di¬ 
agonal matrix of the feedback control gain, and t; is the 
delay time of the ith oscillator’s feedback loop. We as- 
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sume that e > 0 is a small parameter. Each uncoupled 
oscillator x, ; = f, (xj) has the stable limit-cycle solution 
£i{t + Ti) = £i(t) with the natural frequency 17 j = 27t/T- 
We are interested in the case when dissimilarity of the 
periods is of the order of e, i.e., (T — Tj) ~ e for any 
i,j = l,...,N. Also, we assume that the delay times are 
close to the periods, i.e., (T — Tj) ~ e. 

By expanding delayed vector x , : (t — Tj) into a Taylor 
series and omitting higher-than-e-order terms, we arrive 
at the following expression for the control force: 

Kj [Xj(f - Tj) - Xj(£)[ KS 

K i [xj (t - T) - Xj(t)] - KjXj(t - T)AT, (2) 

where AT = Tj — T is the mismatch of the time-delay. 
The first term on the right hand side (r.h.s.) of expres¬ 
sion © is familiar from the controlling chaos, where it 
is used to stabilize unstable periodic orbits in chaotic 
systems [2^, [2^]. Therefore, we use well-known results 
such as the odd number limitation theorem [30| and mis¬ 
matched control scheme pd| . 

The oscillators under DFC, 


x l = f, (xj) + Kj [xj (t - T) - Xj(£)] (3) 


also have the same periodic solutions £j(t) as the free 
oscillators, but with different stability properties, and 
as a consequence, with different perturbation-induced 
phase response. If the free oscillator has an infinites¬ 
imal phase-response curve (iPRC) Zj(i?j) [the iPRC is 
a T-periodic solution of the adjoint equation z ?: (f) = 
— [Dfj(£j(f))] T Zj(t) with the initial condition zf (0) • 
T(0) = 1; see pj-Q] then the oscillator under DFC © 
has the iPRC z| DFC ' l ('i?j) of the same form but with dif¬ 
ferent amplitude (26j: z- DFC ^(dj) = ajZj(dj). The factor 
CKj can be expressed as: 


= aj(Kj) = 


1 + ^ 

m —1 


(4) 


Here the coefficients are the following integrals: 

Cj™) = fj 1 ' z- m \s)t;j m \s)ds, where the upper indices 
denote the particular components of the vectors Zj 
and £j. 

Now we apply the phase reduction technique [26j, [13 
to the oscillator network © assuming that the unper¬ 
turbed oscillators are described by equations © and that 
the perturbation contains two parts: the interoscillatory 
coupling terms G z] and the second term on the r.h.s. of 
expression ©. Both parts are of the same order: O(e). 
The equations for the phase dynamics are 

N 

di = 1 +£Qj(K i )zf(d l )^G,y (£j(0j 
j =1 

-ATaj(Kj)zf(^j)Kj4j(^j). (5) 


Here in the last term of the r.h.s. we write £j (i?j(t)) 
instead of £j($j(i — Tj)). It can be done because 

£i(Mt - Tj)) = + 0(e)) = + 0(e) (6) 

and after multiplication by AT we get the second order 
correction 0(e 2 ) which is omitted [31| . 

The equations for the phase dynamics © are valid only 
when all periodic solutions &(£) are stable solutions of 
the system©. According to the odd number limitation 
theorem [30| , the periodic solution £j(f) is the unstable 
solution of the system © if the condition 


E < -1 (7) 

m —1 


holds. The condition © shows which values of the 
feedback control gains cannot be correctly described by 
Eq. ©. If this condition does not hold, then it is still 
not guaranteed that the solution £i(t) is stable and that 
Eq. © is valid. However, as we see below, for particular 
systems (namely, diffusively-coupled Stuart-Landau and 
FitzHugh -Nagumo models) the condition © is necessary 
and sufficient. 

The phases $j(t) in equation © vary from 0 to T- 
However, when we investigate phase synchronization, it is 
convenient to have phases growing from 0 to 27r. Further¬ 
more, on the r.h.s. of Eq. ©, the first term corresponds 
to trivial growing of the phase. Therefore, we introduce 
new variables = Hj , dj(£) — =^-t, where T is the so- 
called “averaged” period. The number T is not necessar¬ 
ily equal to the average of all oscillator periods and can 
be chosen freely with one requirement: (T — T) ~ £ for 
all i = 1,..., IV. In the new variables, Eq. © can be 
written as 


<Pi=UJi+ £fijO!j(Kj)zf ( ^ ^ 


N 

X E Gi i ( & 

j =i 


n, + a 1 


** + S-t 


j J 


-ATOjaj(Kj)z? 


T ( ''Pi 


n 


Hj Hj 


t Kj& 


Ti + — t 1 

fij a, ' 1 ’ 


where 12 = 2n/T is the “averaged” frequency and u>i = 
fli — 12. The r.h.s. of the last equation depends peri¬ 
odically on time with period T; also, Wj and AT are 
small parameters. Thus we can apply the averaging 
method [321, [33[. Let us denote the averaged phases ipi ( t ). 
The phase model for the averaged phases is 


N 




( 8 ) 


i=i 


where we introduce effective coupling strengths 


eati (Kj) 


(9) 
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effective frequencies 

AT 

uf = Ui + tt-jT [<*» ( K 0 - !], ( 10 ) 


[in equation (flUl) the index i near the U and T is skipped 
without loss of accuracy] and coupling functions 


H ij(X ) = 77T 
-L i 



Hereafter we assume that all network units are 
near identical and described by similar equations, i.e., 
[f» (&(£)) — f ~ e for all*, j = 1... N. We denote 
“averaged” oscillator as x = f (x) , which has stable peri¬ 
odic solution £(i+T) = £(i) and the corresponding iPRC 
z(t + T) = z{t). The choice of the “averaged” oscillator 
must satisfy one requirement: [f(£(f)) — f* (&(£))] ~ £ 
for all i = 1... N. Also, we assume that interoscillatory 
functions have the form G ij (xj, Xj) = a^g (xj, xj) where 
the coefficients a ,; 7 > 0 play the role of the network’s 
adjacency matrix elements. We consider the undirected 
network, therefore the adjacency matrix A T = A. Now 
we can simplify the network’s phase model ([5)1 . Since 
£i(s/fli) = g(s/(l) + 0(e) and z= z(s/fi)+ 0(s), 
without loss of accuracy, the indices i near a can be 
dropped and Eq. ([ 8 ]) can be rewritten as 


ipi = uf + sf ^2 a k 

i=i 


Ki>i ~ V’i). 


( 12 ) 


with effective coupling strengths e® ff = ea(Kj), effective 
frequencies wf® = cj, + (ATj/T) [a(K.j) — 1], coupling 
functions 


Kx) = j; 

r 2n 


(s) 


l 1 

VO 


X- 


o 


ds , (13) 


-l 


and the factor a(K,) = 1 + l , where 

the integrals C= / Q T z^ m \s)^ m \s)ds. 

Our main result is the phase model (fl2l) . Using m, 
we can analyze three important control regimes. Before 
that let us make two additional assumptions: 1 ) - we as¬ 
sume that the condition 0 is necessary and sufficient, 
i.e., the periodic solutions (t) are stable till 0 does 
not hold, and 2 ) - we assume that for uncontrolled net¬ 
work (Kj = 0) there exist a positive threshold coupling 
strength e th > 0 such that when e > £ th the network 
possess stable phase synchronization regime 


ipi=ip2 = ■■■ = 4 >n, 


(14) 


while for e < e t h the phase synchronization cannot be 
achieved. Hence the important control cases are 


(i) If all delay times and feedback control gains are 
equal (r,; = r and Kj = K for all i = 1... N), then syn¬ 
chronization of the network cannot be controlled. In this 
case the phase model m is equivalent to phase model of 
uncontrolled network. It can be seen, if we choose “aver¬ 
aged” period T = r (without loss of generality we always 
can do that) and rewrite effective frequencies 

AT 

w f = Wt + nyKK)-i] 

= u>i + LOi [a(K) — 1] + 0(e 2 ) « WjO:(K). (15) 

Since £® ff = £a(K), the factor a(K) can be eliminated 
from Eqs. m by time-scaling transformation. 

(ii) If all delay times are equal to the periods (r, = Tj 
for all* = 1... N), then = Wj, and it is possible to 
synchronize or desynchronize network independently on 
how big or small the natural interoscillatory coupling is. 
Let us say that all diagonal matrices K, have only one 
nonzero element k\ [> = K ^ and assume that C ^ is 
positive. Then, from condition 0 the phase model (fT 2 l) 
is valid, if K W is in the interval (—l/C^ 1 ), + 00 ). The 
effective coupling strengths £® ff = £/[l + K^C^] go to 
infinity if K W —> —1/C ^ +0. Therefore, the oscillators’ 
phases become extremely sensitive to any perturbations. 
At boundary K ^ = — l/C^ 1 ) all oscillators become neu¬ 
trally stable. For the other boundary, if K W —>• +00 
then £f ff —>• 0 , and the oscillators’ phases cannot syn¬ 
chronize to each other. Note that even if a magnitude of 
the effective coupling strength can be chosen freely, the 
sign cannot be changed. 

(iii) If the mismatch times are equal to 


A Ti LOi 

— ~ fi[l-a(Ki)]’ 


(16) 


then iof s = 0. In this case the network possesses stable 
full phase synchronization regime 


=ip 2 (t) = ...= ipN(t) = const (17) 

under additional conditions: h{ 0 ) = 0 and h'(0) = 7 > 0 . 
The condition h( 0) = 0 is always fulfilled if the units are 
diffusively-coupled, and the condition 7 > 0 represents 
attractive coupling. The stability of the full phase syn¬ 
chronization regime C3 is determined by eigenvalues of a 
matrix M = — 7EL, where E = diag[£f®, £ 3 ®,..., £|f] is a 
diagonal positive-definite matrix and L = D A is a net¬ 
work Laplacian matrix (here D = diag[di, d 2 , ■ ■., d/v] is 
the degree matrix with the elements d, = XOjLi a b')- 
assume that the network is connected and unidirected, 
so the matrix L has eigenvalues 0 = Ai < < ■ ■ ■ < \n- 

By defining a square root of the matrix E as E 1 / 2 
with the entries (ef ®) 1 ^ 2 on the diagonal, we can see 
that M has the same eigenvalues as a symmetric matrix 
M 7 = — 7E 1 / 2 LE 1 / 2 . M' is negative semidefinite matrix 
with only one eigenvalue equal to zero, which corresponds 
to a shift of all phases by a same amount. Hence the full 
phase synchronization regime is stable. 
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FIG. 1. (Color online) Control of synchronization in the 
SL network m for the mismatches AT) = 0. (a) and 

(b) “local” periods, (c) and (d) Kuramoto order parame¬ 
ter, (e) and (f) DFC force applied to the first unit Fi(t) = 
K [x\(t — T\) — xi(t)). The vertical dashed red line shows 
time moment, when the control is turned-on. (a), (c), (e) 
DFC causes synchronization in network with parameters e = 
9 x 10 -4 and K = —0.3, since e eff > e t h- (b), (d), (f) 
DFC causes desynchronization in network with parameters 
e = 5 x IGF 2 and K = 4, since £ eS < eth- 


III. NUMERICAL DEMONSTRATIONS 

As a first example, we study a network of N = 8 all- 
to-all diffusively coupled Stuart- Landau (SL) models de¬ 
scribed by the following equations: 


N 

±i = Xi (l — x? — Ui) - SliUi + eN~ l ^2 2 (xj - Xi) 

0 =1 

+K [Xi(t - Ti) - Xi(t)] , (18a) 

Vi = Vi (l -Vi) +ViXi- (18b) 

As an “averaged” oscillator we choose SL model with 
0 = 1. For this case, the periodic orbit and iPRC 
can be found analytically: £(f) = [cost, sin t] T 
and z(i) = [—sint,cost] T . According to (fT3l) . 

the coupling function h{\) = sin y. For the 

oscillators natural frequencies we choose follow¬ 
ing values Ot = O + u>i, where Wj = 10~ 3 x 
{1.38, 2.54, -1.93, -4.87, -2.12, 3.95, 4.31, -3.26}. 

For SL model the factor a(K) = [1 + Kt:]^ 1 and 
the feedback control gain K can be selected from the 
interval (—7r -1 , +oo). We calculate numerically that 
the network m without control (K = 0) possesses 
phase synchronization if e is above the threshold value 
eth = 7x KT 3 . The adjacency matrix of the network (fTSl) 
is dij = N _1 and the corresponding phase model is the 
celebrated Kuramoto model 

N 

tpi = + £ eff /V -1 ^2 sin(i pj - ipi). (19) 

j =i 
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FIG. 2. (Color online) The synchrony-asynchrony transition 
in the SL network for the control regime A Ti = 0. A snap¬ 
shot of the “local” periods versus the control gain calculated 
numerically for different coupling strengths: (a) e = 9 x 1CU 4 
and (b) e = 5 x 10~ 2 . A vertical red line shows analytically 
derived critical value of the control gain. 


As a synchronization criteria we choose two mea¬ 
surements: Kuramoto order parameter r[t) = 

IV -1 Y^j=i ex Pand the “local” periods Xi oca i (or 


sometimes called interspike intervals) defined as a time 
interval between two neighboring maxima of the first dy¬ 
namical variable. Figure [T] shows numerical simulation 
of the SL network m when the mismatches AT) = 0. 
The transition synchrony-asynchrony occurs at a criti¬ 
cal control gain K c = 7 r -1 [e/e t h — 1]. In Figure. [2] we 
demonstrate the synchrony-asynchrony transition in the 
SL network. As we can see, analytical results coincide 
with numerical simulations. 


In order to demonstrate full phase synchronization 
regime, we simulate SL network with the same param¬ 
eters as presented in Figs. Q] (a), (c) and (e), only the 
mismatch times are selected according to m- Th e “lo¬ 
cal” periods and F\ are very similar to that presented in 
Figs. [I] (a) and (e), only the Kuramoto order parameter 
is much closer to one (cf. Fig. [3]). 



time 

FIG. 3. (Color online) Kuramoto order parameters of SL 
network (1181) depicted on a semilog plot. The blue thin line 
reproduces results from Fig. |T] (c), while the green thick line 
calculated in full-phase-synchronization regime. The verti¬ 
cal dashed red line shows time moment, when the control is 
turned-on. 


In order to show that analytical results valid for non¬ 
trivial oscillators and for the nontrivial network topology, 
we investigate the network of N = 8 diffusivclycoupled 
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FIG. 4. Topology of the FitzHugh-Nagumo oscillator net¬ 
work m - 



FIG. 5. (Color online) Control of synchronization in the FHN 
network 1201) for the mismatches AT) = 0. (a) and (b) “local” 
periods, (c) and (d) DFC force applied to the first unit F\ (t) = 
K [xi(t — Ti) — xi(t)]. The vertical dashed red line shows 
time moment, when the control is turned-on. (a), (c) DFC 
causes synchronization in network with parameters e = 5 x 
10 -5 and K = —0.09, since e eff > eth- (b), (d) DFC causes 
desynchronization in network with parameters e — 10~ 3 and 
K = 0.5, since e eff < eth- 


FitzHugh -Nagumo (FHN) [34], HH models 

N 

Xi = Xi - x?/3 — yi + 0.5 + e'^ 2 a ij( x j - x i ) 

3 =1 

+K [Xi (t - Ti ) - Xi (t )] , (20a) 

Vi = £i(xi + 0.7 — 0.8j/i). (20b) 

The adjacency matrix elements a,ij = 1, if the unit i is 
connected to the unit j, and a^- = 0 otherwise. The 
network topology is illustrated in Fig. [4] As an “aver¬ 


aged” oscillator we choose FHN model with e = 0.08. 
For such model the constant C W ss 10.02 computed nu¬ 
merically. We check that “averaged” oscillator possesses 
stable periodic solution £(i), when control gain K is in 
the interval (—1/©© +oo). In the network (l20l) each os¬ 
cillator has different parameter e, = e +Ae*, where Ae^ = 
10” 1 * 3 4 5 6 7 x {0.3, -1.7, -0.9, 2.1, 1.5, -2.6, -1.1, 0.8}, and 
without control (K = 0) it possesses phase synchroniza¬ 
tion if £ is above the threshold £ t h = 3.6 x 10 -4 . Figure[5] 
shows numerical simulation of the FHN network EU1) 
when the mismatches AT) = 0. Again, analytical results 
coincide with numerical simulations. 


IV. CONCLUSION 

We present framework for controlling synchrony in 
weakly coupled oscillator networks by delayed feedback 
control. We show that when the delay time is close 
to the period of a particular oscillator, the network’s 
phase model almost coincides with the uncontrolled net¬ 
work’s phase model. The only difference is that effective 
coupling strengths and effective frequencies depend on 
control parameters. By appropriate choice of the con¬ 
trol parameters the magnitude of the effective coupling 
strength can be selected arbitrary, while the sign cannot 
be changed. Unlike coupling strength, the sign of the 
effective frequencies can be inverted. 

In this work we have restricted ourselves to the case 
when control term appears as an external force applied 
to the oscillator. However it can be simply generalized 
to the case of arbitrary functional dependence of the os¬ 
cillator on the control signal. 
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